	! Unit cell utilities

	module cellutils
	implicit none
	contains

	subroutine fold(cell,icell,x1,y1,z1)
	implicit none
	real*8, intent(in) :: cell(9),icell(9)
	real*8, intent(inout) :: x1,y1,z1
	! Multiply coordinate by inverse cell to acquire unit cell coordinates
	call vmatmult(x1,y1,z1,icell)
	x1 = x1 - NINT(x1)
	y1 = y1 - NINT(y1)
	z1 = z1 - NINT(z1)
	call vmatmult(x1,y1,z1,cell)
	end subroutine fold

	subroutine mimvec(cell,icell,x1,y1,z1,x2,y2,z2,x3,y3,z3)
	implicit none
	! Performs minimum image convention.
	! xyz1: point to consider; xyz2: reference point; xyz3: minimum image coords of xyz1 (result)
	real*8, intent(in) :: x1,y1,z1,x2,y2,z2,cell(9),icell(9)
	real*8, intent(out) :: x3,y3,z3
	real*8 :: x,y,z
	! Multiply coordinate by inverse cell to acquire unit cell coordinates
	x = x1 - x2
	y = y1 - y2
	z = z1 - z2
	call vmatmult(x,y,z,icell)
	x = x - NINT(x)
	y = y - NINT(y)
	z = z - NINT(z)
	call vmatmult(x,y,z,cell)
	x3 = x
	y3 = y
	z3 = z
! 	x3 = x1 - x2
! 	y3 = y1 - y2
! 	z3 = z1 - z2
	end subroutine mimvec

	subroutine calc_icell(icell)
	use inverse
	implicit none
	real*8, intent(inout) :: icell(9)
	call gj(icell,3,3)
	end subroutine calc_icell

	subroutine getvector(x1,y1,z1,x2,y2,z2,x3,y3,z3)
	implicit none
	real*8 :: x1,y1,z1,x2,y2,z2,x3,y3,z3
	! xyz1,xyz2: coordinates of two atoms; xyz3: vector between them (result)
	x3=x1-x2
	y3=y1-y2
	z3=z1-z2
	end subroutine getvector

	subroutine vmatmult(x1,y1,z1,a)
	implicit none
	real*8 :: x1,y1,z1,x2,y2,z2,a(9)
	! Multiply vector(3) v by matrix(3x3) a
	! Assume matrix is stored in row-major order
	x2 = x1*a(1) + y1*a(4) + z1*a(7)
	y2 = x1*a(2) + y1*a(5) + z1*a(8)
	z2 = x1*a(3) + y1*a(6) + z1*a(9)
	x1 = x2
	y1 = y2
	z1 = z2
	end subroutine vmatmult

	real*8 function volume(cell)
	implicit none
	real*8, intent(in) :: cell(9)
	! Calculate determinant of axes matrix
	volume = 0.0d0
	volume = volume + cell(1) * (cell(5)*cell(9) - cell(8)*cell(6))
	volume = volume - cell(2) * (cell(4)*cell(9) - cell(7)*cell(6))
	volume = volume + cell(3) * (cell(4)*cell(8) - cell(7)*cell(5))
	end function volume

	end module cellutils
